Data Set-up

#############################################################################################
setwd("~/Dropbox/CalcanealGearRatio")
calcmeasure <- read.csv("~/Dropbox/CalcanealGearRatio/LACM_CalcaneusMeasurements_2018_8_29.csv", stringsAsFactors=FALSE)
calcmeasure <- calcmeasure[calcmeasure["Measurement.Type.Manual.Photo."] == "Photo",1:14]
calcSpecies <- unique(paste(calcmeasure$Genus, calcmeasure$Species,sep="_"))

MOM.Mat <- read.csv("~/Dropbox/CalcanealGearRatio/MOM_v4.1.csv")
#remove empty rows
MOM.Mat <- MOM.Mat[!MOM.Mat$Species == "",]
#aggregate to get single species entries
MOM.Mat[,"TaxonomicName"] <- paste(MOM.Mat[,"Genus"],MOM.Mat[,"Species"], sep="_")
#MOM.Mat[,"TaxonomicName"] <- gsub(pattern = " ", replacement =  "_", x = MOM.Mat[,"TaxonomicName"],)

MOM.Mat[,"CombinedMass(kg)"] <- MOM.Mat[,"Combined.Mass..g."]/1000
MOM.Mat[,"logMass.Kg"] <- log10(MOM.Mat[,"CombinedMass(kg)"])
MOM.Mat <- MOM.Mat[,c("Order","FAMILY","Genus","Species","TaxonomicName","CombinedMass(kg)","logMass.Kg")]
#combine same species if multiple entires(e.g. Muntiacus_reevesi has separate entries for mainland Asia and Insular)
MOM.Mat <- aggregate(MOM.Mat[,-c(1:5)], by = list(MOM.Mat$Order, MOM.Mat$FAMILY,MOM.Mat$Genus,
                                                  MOM.Mat$Species, MOM.Mat$TaxonomicName), mean)
colnames(MOM.Mat)[1:5] <- c("Order","FAMILY","Genus","Species","TaxonomicName")

ecoMat <- read.csv("~/Dropbox/CalcanealGearRatio/ExtantUngulate_BM_Eco_Mendoza_etal2007.csv")
ecoMat["TaxonomicName"] <- paste(ecoMat[,"Genus"],paste("_",ecoMat[,"Species"],sep=""),sep="")

#############################################################################################
###Analysis
#############################################################################################
getOneSpecimenAverage <- function(this.specimen, calcmeasure) {
    this.measurements <- c("Length.Calc", "Length..medial.", "Length..lateral.", "Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.", "Length.to.Sustentacular.Process", "Length.to.Calcaneoastragular.facet", "Calcaneoastragular.facet.to.cuboid.facet")
    colMeans(calcmeasure[calcmeasure$Catalog.Number== this.specimen, this.measurements], na.rm=TRUE)
}

getSpecimenAverages <- function(calcmeasure) {
    specimen.vec <- sort(unique(calcmeasure$Catalog.Number))
    n <- t(sapply(specimen.vec, function(x) calcmeasure[calcmeasure$Catalog.Number==x, c("Museum.Acronym", "Catalog.Number", "Genus", "Species")][1,]))
    m <- data.frame(n, t(sapply(specimen.vec, getOneSpecimenAverage, calcmeasure)))
    m
}

getSpeciesAverages <- function(measureMat) {
  rownames(measureMat) <- NULL
  species.vec <- sort(unique(measureMat$taxon))
  this.measurements <- c("Length.Calc", "Length..medial.", "Length..lateral.", "Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.", "Length.to.Sustentacular.Process", "Length.to.Calcaneoastragular.facet", "Calcaneoastragular.facet.to.cuboid.facet")
  tst <- species.vec[1]
  this.species <- tst
  spMat <-matrix(nrow=length(species.vec),ncol = length(this.measurements))
  for(jj in seq(1, length(species.vec),1))
  {
    this.species <- species.vec[jj]
    spMat[jj,] <- colMeans(measureMat[measureMat$taxon == this.species, this.measurements], na.rm=TRUE)
  }
  colnames(spMat) <- this.measurements
  rownames(spMat) <- species.vec
  spMat
}

m <- getSpecimenAverages(calcmeasure)
m$taxon <- paste(m$Genus, tolower(m$Species), sep="_")
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Calcaneoastragular.facet.to.cuboid.facet #this is based on the figure caption
m$R2 <- m$Calcaneoastragular.facet.to.cuboid.facet/m$Length.to.Calcaneoastragular.facet
#m$R2  <- m$Length.to.Calcaneoastragular.facet/m$Length.Calc #this is based on the designation listed in the figure
m$R3 <- m$Length.Calc/m$Length.to.Sustentacular.Process
m$A1 <- rowMeans(m[,c("Length..medial.", "Length..lateral.")], na.rm=TRUE)/rowMeans(m[,c("Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.")], na.rm=TRUE)
m$habitat <- ecoMat$Habitat[match(m$taxon, ecoMat$TaxonomicName)]
m$BM <- MOM.Mat$logMass.Kg[match(m$taxon, MOM.Mat$TaxonomicName)]

#family
m$Family <- MOM.Mat$FAMILY[match(m$taxon, MOM.Mat$TaxonomicName)]
#order
#family
m$Order <- MOM.Mat$Order[match(m$taxon, MOM.Mat$TaxonomicName)]

##remove Perissodactyls (3 species)
m <- m[!m$Order == "Perissodactyla",]

##species level
m.sp <- as.data.frame(getSpeciesAverages(m))
m.sp$taxon <- row.names(m.sp)
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Calcaneoastragular.facet.to.cuboid.facet #this is based on the figure caption
m.sp$R2 <- m.sp$Calcaneoastragular.facet.to.cuboid.facet/m.sp$Length.to.Calcaneoastragular.facet
#m$R2  <- m$Length.to.Calcaneoastragular.facet/m$Length.Calc #this is based on the designation listed in the figure
m.sp$R3 <- m.sp$Length.Calc/m.sp$Length.to.Sustentacular.Process
m.sp$A1 <- rowMeans(m.sp[,c("Length..medial.", "Length..lateral.")], na.rm=TRUE)/rowMeans(m.sp[,c("Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.")], na.rm=TRUE)
m.sp$habitat <- ecoMat$Habitat[match(m.sp$taxon, ecoMat$TaxonomicName)]
m.sp$BM <- MOM.Mat$logMass.Kg[match(m.sp$taxon, MOM.Mat$TaxonomicName)]

#reorder
m.sp <- m.sp[,c(9,13,1:8,14, 10:12)]
m.sp <- m.sp[complete.cases(m.sp),] #68 species
m <- m.sp
#############################################################################################
habitat.colors <- c("dodgerblue", "chartreuse", "goldenrod1")

Polly 2010 Figure 4c

    # habitat.colors <- rainbow(length(unique(m$habitat)))
    plot(m$R2, m$R3, xlim=c(0.35,0.65), ylim=c(1.1, 1.5), pch=21, cex=1, bg=habitat.colors[m$habitat], xlab = "Gear Ratio (R2)", ylab = "Gear Ratio (R3)")
    # plot(m$R2, m$R3, pch=21, cex=1, bg=habitat.colors[m$habitat])
    text(m$R2, m$R3, labels=m$taxon, pos=3, cex=0.3)
    abline(lm(m$R3 ~ m$R2))

Polly 2010 Figure 6

    plot(sort(m$R3), pch=21, bg=habitat.colors[m$habitat], type="n", ylim=c(1,1.5))
    rect(xleft=-10, ybottom=mean(m$R3, na.rm=TRUE)-sd(m$R3, na.rm=TRUE), xright=150, ytop=mean(m$R3, na.rm=TRUE)+sd(m$R3, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
    abline(h=mean(m$R3, na.rm=TRUE), lty=2)
    points(sort(m$R3), pch=21, bg=habitat.colors[m$habitat[order(m$R3)]])

    #text(seq_len(nrow(m)), sort(m$R3), labels=m$taxon[order(m$R3)], pos=3, cex=0.2)
    #axis(1, at = seq_len(nrow(m)), labels = sort(m$R3)) make diagonal labels for each taxon

Polly 2010 Figure 6 but for R2, instead of R3

    plot(sort(m$R2), pch=21, bg=habitat.colors[m$habitat], type="n", ylab = "Gear Ratio (R2)")
    rect(xleft=-10, ybottom=mean(m$R2, na.rm=TRUE)-sd(m$R2, na.rm=TRUE), xright=150, ytop=mean(m$R2, na.rm=TRUE)+sd(m$R2, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
    abline(h=mean(m$R2, na.rm=TRUE), lty=2)
    points(sort(m$R2), pch=21, bg=habitat.colors[m$habitat[order(m$R2)]])
    text(seq_len(nrow(m)), sort(m$R2), labels=m$taxon[order(m$R2)], pos=3, cex=0.3)

Polly 2010 Figure 6 but for A1, instead of R3

    plot(sort(m$A1), pch=21, bg=habitat.colors[m$habitat], type="n", ylab = "Area of Astragalus")
    rect(xleft=-10, ybottom=mean(m$A1, na.rm=TRUE)-sd(m$A1, na.rm=TRUE), xright=150, ytop=mean(m$A1, na.rm=TRUE)+sd(m$A1, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
    abline(h=mean(m$A1, na.rm=TRUE), lty=2)
    points(sort(m$A1), pch=21, bg=habitat.colors[m$habitat[order(m$A1)]])
    text(seq_len(nrow(m)), sort(m$A1), labels=m$taxon[order(m$A1)], pos=3, cex=0.3)

### Correlation of ratios with body mass

    this.ratio <- m$R2
    plot(m$BM, this.ratio, pch=21, bg=habitat.colors[m$habitat], ylab = "Gear Ratio (R2)", xlab = "Body mass log(kg)")
    text(m$BM, this.ratio, labels=m$taxon, pos=3, cex=0.3)
    abline(lm(this.ratio ~ m$BM))

    cor.test(m$BM, this.ratio)
## 
##  Pearson's product-moment correlation
## 
## data:  m$BM and this.ratio
## t = -0.47015, df = 66, p-value = 0.6398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2921755  0.1831732
## sample estimates:
##        cor 
## -0.0577751
    this.ratio <- m$R3
    plot(m$BM, this.ratio, pch=21, bg=habitat.colors[m$habitat], ylab = "Gear Ratio (R3)", xlab = "Body mass log(kg)")
    text(m$BM, this.ratio, labels=m$taxon, pos=3, cex=0.3)
    abline(lm(this.ratio ~ m$BM))

    cor.test(m$BM, this.ratio)
## 
##  Pearson's product-moment correlation
## 
## data:  m$BM and this.ratio
## t = -1.2026, df = 66, p-value = 0.2334
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.37187876  0.09531537
## sample estimates:
##        cor 
## -0.1464376

3-Habitat Hypothesis

PCA

library(ggfortify)
## Loading required package: ggplot2
dropCol <- c(1:3,11:14)
df <- log(m[,-dropCol])
habitat <- m$habitat

df <- scale(df)

df.pca <- prcomp(df)
summary(df.pca)
## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.6306 0.20277 0.13787 0.10559 0.07592 0.04685
## Proportion of Variance 0.9886 0.00587 0.00272 0.00159 0.00082 0.00031
## Cumulative Proportion  0.9886 0.99442 0.99714 0.99873 0.99955 0.99987
##                            PC7
## Standard deviation     0.03068
## Proportion of Variance 0.00013
## Cumulative Proportion  1.00000
#plot(df.pca$x[,1], df.pca$x[,2], pch = 16, col = habitat.colors[m$habitat], ylab = "PC2", xlab = "PC1")

autoplot(df.pca, data = m, colour = 'habitat', label = FALSE, label.size = 3, 
         loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3)+
  scale_color_manual(breaks = c("CH", "MH", "OH"),
                     values=c("dodgerblue", "chartreuse", "goldenrod1"))

#plot(df.pca$x[,2], df.pca$x[,3], pch = 16, col = habitat.colors[m$habitat], ylab = "PC3", xlab = "PC2")

autoplot(df.pca, data = m, x=2, y=3, choices = c(2:3), colour = 'habitat', 
               label = TRUE, label.size = 3, 
             loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3) +
  scale_color_manual(breaks = c("CH", "MH", "OH"),
                     values=c("dodgerblue", "chartreuse", "goldenrod1"))

K-means

set.seed(1)
#set-symbols for pch
ksym <- gsub("OH",16, m$habitat)
ksym <- gsub("MH",17, ksym)
ksym <- gsub("CH",18, ksym)
ksym <- as.numeric(ksym)

df.kmean <- kmeans(df, 3)

kcol <- gsub(1, "red", df.kmean$cluster)
kcol <- gsub(2, "darkgreen", kcol)
kcol<- gsub(3, "blue", kcol)

autoplot(df.kmean, data = m, size = 0.1) +
  geom_point(shape = ksym, colour = kcol, size = 2) 

Discriminant Analysis

Linear Discriminant Analysis

library(lfda)
library(MASS)

lda.tst <- lda(x = as.matrix(df), grouping = m$habitat)
plot(lda.tst, col = habitat.colors[m$habitat])

Local Fisher Discriminant Analysis (LFDA)

model1_3h <- lfda(df, m$habitat, 3, metric="plain")
autoplot(model1_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
  geom_point(shape = ksym, size = 2) 

autoplot(model1_3h, data = m, x=2, y=3,frame = TRUE, frame.colour = 'habitat', size=0.1)+
  geom_point(shape = ksym, size=2) 

###Kernal Local Fisher Discriminant Analysis (KLFDA)

model2_3h <- klfda(kmatrixGauss(df), m$habitat, 3, metric="plain")
autoplot(model2_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
  geom_point(shape = ksym, size = 2) 

autoplot(model2_3h, data = m,x=2,y=3, frame = TRUE, frame.colour = 'habitat', size=0.1)+
  geom_point(shape = ksym, size = 2) 

Semi-supervised Local Fisher Discriminant Analysis (SELF)

model3_3h <- self(df, m$habitat, beta = 0.1, r = 3, metric="plain")
autoplot(model3_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
  geom_point(shape = ksym, size = 2) 

require(phytools)
## Loading required package: phytools
## Loading required package: ape
## Loading required package: maps
#read in tree
testTree <- read.nexus("/Users/emdoughty/Downloads/Appendix 3 alignment.nex")
testTreeTips <- testTree$tip.label

getTipTaxon <- function(treetips){
  tipList <- strsplit(treetips,split = "_")
  tipMat <- matrix(nrow = length(tipList), ncol = 2)
  for(ii in seq(1, length(tipList),1))
  {
    tipMat[ii,1] <- tipList[[ii]][1]
    tipMat[ii,2] <- tipList[[ii]][2]
  }
  
  labelTipsNew <- paste(tipMat[,1], tipMat[,2], sep="_")
  return (labelTipsNew)
}
labelTipsNew <- getTipTaxon(testTree$tip.label)
#how many entries are present per species name
new.table <- table(labelTipsNew)
removeVec <- vector()

for(yy in seq(1, length(new.table),1))
{
  if(new.table[yy] > 1)
  {
    names(new.table[yy])
    tip.options <- which(labelTipsNew == names(new.table[yy]))
    tip.remove <- sample(tip.options,size = length(tip.options)-1)
    removeVec <- append(removeVec, tip.remove)
  }
}

#check removeVec against names
testTree$tip.label[removeVec]
##  [1] "Alcelaphus_buselaphus_buselaphus"          
##  [2] "Balaenoptera_edeni_brydei_NC006928"        
##  [3] "Bos_javanicus_javanicus"                   
##  [4] "Bos_taurus_taurus_EU177832"                
##  [5] "Bubalus_bubalis_bubalis_AF547270"          
##  [6] "Cephalophus_callipygus_1"                  
##  [7] "Cephalophus_callipygus_2"                  
##  [8] "Cervus_nippon_taiouanus_NC008462"          
##  [9] "Cervus_nippon_centralis_NC006993"          
## [10] "Connochaetes_taurinus_2"                   
## [11] "Dama_dama_dama"                            
## [12] "Eudorcas_rufifrons_2"                      
## [13] "Gazella_dorcas_pelzelnii"                  
## [14] "Gazella_gazella_erlangeri"                 
## [15] "Gazella_subgutturosa_marica"               
## [16] "Giraffa_camelopardalis_angolensis_NC012100"
## [17] "Mazama_americana_2"                        
## [18] "Mazama_nemorivaga_1"                       
## [19] "Naemorhedus_griseus_griseus_FJ207532"      
## [20] "Odocoileus_virginianus_3"                  
## [21] "Odocoileus_virginianus_1"                  
## [22] "Pecari_tajacu_2_NC012103"                  
## [23] "Pecari_tajacu_3"                           
## [24] "Philantomba_monticola_2"                   
## [25] "Tragelaphus_scriptus_2"                    
## [26] "Tragelaphus_scriptus_3"
#collapse/remove species with multiple entries
testTree_crop <- drop.tip(testTree, removeVec)

#check if it worked
newNameList <- getTipTaxon(testTree_crop$tip.label)
table(newNameList)
## newNameList
##        Addax_nasomaculatus         Aepyceros_melampus 
##                          1                          1 
##      Alcelaphus_buselaphus                Alces_alces 
##                          1                          1 
##          Ammotragus_lervia     Antidorcas_marsupialis 
##                          1                          1 
##      Antilocapra_americana        Antilope_cervicapra 
##                          1                          1 
##       Arabitragus_jayakari                  Axis_axis 
##                          1                          1 
##              Axis_porcinus         Balaena_mysticetus 
##                          1                          1 
## Balaenoptera_acutorostrata   Balaenoptera_bonaerensis 
##                          1                          1 
##      Balaenoptera_borealis         Balaenoptera_edeni 
##                          1                          1 
##      Balaenoptera_musculus        Balaenoptera_omurai 
##                          1                          1 
##      Balaenoptera_physalus          Berardius_bairdii 
##                          1                          1 
##                Bison_bison              Bison_bonasus 
##                          1                          1 
##     Blastocerus_dichotomus                 Bos_gaurus 
##                          1                          1 
##              Bos_grunniens              Bos_javanicus 
##                          1                          1 
##                 Bos_taurus    Boselaphus_tragocamelus 
##                          1                          1 
##            Bubalus_bubalis     Bubalus_depressicornis 
##                          1                          1 
##         Budorcas_taxicolor         Camelus_bactrianus 
##                          1                          1 
##        Camelus_dromedarius              Camelus_ferus 
##                          1                          1 
##          Caperea_marginata            Capra_caucasica 
##                          1                          1 
##            Capra_falconeri               Capra_hircus 
##                          1                          1 
##                 Capra_ibex              Capra_nubiana 
##                          1                          1 
##            Capra_pyrenaica             Capra_sibirica 
##                          1                          1 
##        Capreolus_capreolus        Capricornis_crispus 
##                          1                          1 
## Capricornis_milneedwardsii       Capricornis_swinhoei 
##                          1                          1 
##         Cephalophus_adersi     Cephalophus_callipygus 
##                          1                          1 
##       Cephalophus_dorsalis       Cephalophus_jentinki 
##                          1                          1 
##    Cephalophus_leucogaster     Cephalophus_natalensis 
##                          1                          1 
##     Cephalophus_nigrifrons        Cephalophus_ogilbyi 
##                          1                          1 
##      Cephalophus_rufilatus    Cephalophus_silvicultor 
##                          1                          1 
##         Cephalophus_spadix Cephalorhynchus_heavisidii 
##                          1                          1 
##             Cervus_elaphus              Cervus_nippon 
##                          1                          1 
##     Choeropsis_liberiensis          Connochaetes_gnou 
##                          1                          1 
##      Connochaetes_taurinus                  Dama_dama 
##                          1                          1 
##        Damaliscus_pygargus         Delphinus_capensis 
##                          1                          1 
##      Dorcatragus_megalotis      Elaphodus_cephalophus 
##                          1                          1 
##       Elaphurus_davidianus               Equus_asinus 
##                          1                          1 
##      Eschrichtius_robustus        Eubalaena_australis 
##                          1                          1 
##         Eubalaena_japonica         Eudorcas_rufifrons 
##                          1                          1 
##          Gazella_bennettii            Gazella_cuvieri 
##                          1                          1 
##             Gazella_dorcas            Gazella_gazella 
##                          1                          1 
##         Gazella_leptoceros             Gazella_spekei 
##                          1                          1 
##       Gazella_subgutturosa     Giraffa_camelopardalis 
##                          1                          1 
##            Grampus_griseus      Hemitragus_jemlahicus 
##                          1                          1 
##    Hippocamelus_antisensis     Hippopotamus_amphibius 
##                          1                          1 
##        Hippotragus_equinus          Hippotragus_niger 
##                          1                          1 
##         Hydropotes_inermis       Hyemoschus_aquaticus 
##                          1                          1 
##      Hyperoodon_ampullatus           Inia_geoffrensis 
##                          1                          1 
##       Kobus_ellipsiprymnus                Kobus_leche 
##                          1                          1 
##            Kogia_breviceps Lagenorhynchus_albirostris 
##                          1                          1 
##              Lama_guanicoe         Lipotes_vexillifer 
##                          1                          1 
##        Litocranius_walleri             Madoqua_kirkii 
##                          1                          1 
##           Madoqua_saltiana           Mazama_americana 
##                          1                          1 
##         Mazama_gouazoubira          Mazama_nemorivaga 
##                          1                          1 
##              Mazama_rufina     Megaptera_novaeangliae 
##                          1                          1 
##          Monodon_monoceros        Moschus_berezovskii 
##                          1                          1 
##        Moschus_moschiferus       Muntiacus_crinifrons 
##                          1                          1 
##          Muntiacus_muntjak          Muntiacus_reevesi 
##                          1                          1 
##     Muntiacus_vuquangensis        Naemorhedus_baileyi 
##                          1                          1 
##        Naemorhedus_griseus                Nanger_dama 
##                          1                          1 
##              Nanger_granti       Nanger_soemmerringii 
##                          1                          1 
##           Neotragus_batesi       Nesotragus_moschatus 
##                          1                          1 
##        Odocoileus_hemionus     Odocoileus_virginianus 
##                          1                          1 
##           Okapia_johnstoni        Oreamnos_americanus 
##                          1                          1 
##      Oreotragus_oreotragus                 Oryx_beisa 
##                          1                          1 
##                Oryx_dammah               Oryx_gazella 
##                          1                          1 
##              Oryx_leucoryx             Ourebia_ourebi 
##                          1                          1 
##           Ovibos_moschatus                 Ovis_aries 
##                          1                          1 
##    Ozotoceros_bezoarcticus            Panthera_tigris 
##                          1                          1 
##       Pantholops_hodgsonii              Pecari_tajacu 
##                          1                          1 
##            Pelea_capreolus     Phacochoerus_africanus 
##                          1                          1 
##       Philantomba_maxwelli      Philantomba_monticola 
##                          1                          1 
##          Phocoena_phocoena     Physeter_macrocephalus 
##                          1                          1 
##       Platanista_gangetica     Pontoporia_blainvillei 
##                          1                          1 
##       Potamochoerus_porcus         Procapra_gutturosa 
##                          1                          1 
##   Przewalskium_albirostris            Pseudois_nayaur 
##                          1                          1 
##    Pseudoryx_nghetinhensis        Pteropus_dasymallus 
##                          1                          1 
##        Pudu_mephistophiles                  Pudu_puda 
##                          1                          1 
##       Rahicerus_campestris          Rangifer_tarandus 
##                          1                          1 
##          Redunca_arundinum        Redunca_fulvorufula 
##                          1                          1 
##         Rucervus_duvauceli              Rucervus_eldi 
##                          1                          1 
##        Rupicapra_pyrenaica        Rupicapra_rupicapra 
##                          1                          1 
##               Rusa_alfredi            Rusa_timorensis 
##                          1                          1 
##              Rusa_unicolor             Saiga_tatarica 
##                          1                          1 
##            Sousa_chinensis         Stenella_attenuata 
##                          1                          1 
##      Stenella_coeruleoalba               Sus_barbatus 
##                          1                          1 
##                 Sus_scrofa         Sylvicapra_grimmia 
##                          1                          1 
##            Syncerus_caffer    Tetracerus_quadricornis 
##                          1                          1 
##        Tragelaphus_angasii      Tragelaphus_derbianus 
##                          1                          1 
##      Tragelaphus_eurycerus       Tragelaphus_imberbis 
##                          1                          1 
##           Tragelaphus_oryx       Tragelaphus_scriptus 
##                          1                          1 
##         Tragelaphus_spekii   Tragelaphus_strepsiceros 
##                          1                          1 
##           Tragulus_kanchil           Tursiops_aduncus 
##                          1                          1 
##         Tursiops_truncatus              Vicugna_pacos 
##                          1                          1
testTree_crop$tip.label <- newNameList

#find which species are not on the tree
unique(m$taxon[!m$taxon %in% newNameList])
##  [1] "Beatragus_hunteri"        "Cephalophus_zebra"       
##  [3] "Cervus_canadensis"        "Damaliscus_lunatus"      
##  [5] "Hemitragus_hylocrius"     "Kobus_kob"               
##  [7] "Lama_pacos"               "Madoqua_guentheri"       
##  [9] "Neotragus_moschatus"      "Ovis_canadensis"         
## [11] "Ovis_dalli"               "Ovis_orientalis"         
## [13] "Phacochoerus_aethiopicus" "Philantomba_maxwellii"   
## [15] "Raphicerus_campestris"    "Rucervus_eldii"          
## [17] "Taurotragus_oryx"         "Tayassu_pecari"          
## [19] "Tragulus_javanicus"       "Tragulus_napu"           
## [21] "Vicugna_vicugna"
unique(testTree_crop$tip.label[!testTree_crop$tip.label %in% m$taxon])
##   [1] "Pteropus_dasymallus"        "Equus_asinus"              
##   [3] "Panthera_tigris"            "Phacochoerus_africanus"    
##   [5] "Potamochoerus_porcus"       "Sus_barbatus"              
##   [7] "Lama_guanicoe"              "Vicugna_pacos"             
##   [9] "Camelus_dromedarius"        "Camelus_ferus"             
##  [11] "Choeropsis_liberiensis"     "Balaena_mysticetus"        
##  [13] "Eubalaena_australis"        "Eubalaena_japonica"        
##  [15] "Caperea_marginata"          "Balaenoptera_acutorostrata"
##  [17] "Balaenoptera_bonaerensis"   "Eschrichtius_robustus"     
##  [19] "Megaptera_novaeangliae"     "Balaenoptera_physalus"     
##  [21] "Balaenoptera_musculus"      "Balaenoptera_omurai"       
##  [23] "Balaenoptera_borealis"      "Balaenoptera_edeni"        
##  [25] "Kogia_breviceps"            "Physeter_macrocephalus"    
##  [27] "Platanista_gangetica"       "Hyperoodon_ampullatus"     
##  [29] "Berardius_bairdii"          "Lipotes_vexillifer"        
##  [31] "Inia_geoffrensis"           "Pontoporia_blainvillei"    
##  [33] "Monodon_monoceros"          "Phocoena_phocoena"         
##  [35] "Lagenorhynchus_albirostris" "Cephalorhynchus_heavisidii"
##  [37] "Grampus_griseus"            "Stenella_attenuata"        
##  [39] "Sousa_chinensis"            "Tursiops_truncatus"        
##  [41] "Stenella_coeruleoalba"      "Delphinus_capensis"        
##  [43] "Tursiops_aduncus"           "Tragulus_kanchil"          
##  [45] "Okapia_johnstoni"           "Capreolus_capreolus"       
##  [47] "Ozotoceros_bezoarcticus"    "Mazama_gouazoubira"        
##  [49] "Hippocamelus_antisensis"    "Blastocerus_dichotomus"    
##  [51] "Mazama_nemorivaga"          "Mazama_rufina"             
##  [53] "Pudu_mephistophiles"        "Elaphodus_cephalophus"     
##  [55] "Muntiacus_muntjak"          "Muntiacus_crinifrons"      
##  [57] "Muntiacus_vuquangensis"     "Rucervus_duvauceli"        
##  [59] "Elaphurus_davidianus"       "Rucervus_eldi"             
##  [61] "Rusa_alfredi"               "Cervus_elaphus"            
##  [63] "Przewalskium_albirostris"   "Moschus_berezovskii"       
##  [65] "Moschus_moschiferus"        "Tetracerus_quadricornis"   
##  [67] "Tragelaphus_scriptus"       "Tragelaphus_angasii"       
##  [69] "Tragelaphus_oryx"           "Tragelaphus_derbianus"     
##  [71] "Syncerus_caffer"            "Bubalus_bubalis"           
##  [73] "Bubalus_depressicornis"     "Pseudoryx_nghetinhensis"   
##  [75] "Bison_bonasus"              "Bos_taurus"                
##  [77] "Bos_grunniens"              "Bos_gaurus"                
##  [79] "Neotragus_batesi"           "Nesotragus_moschatus"      
##  [81] "Procapra_gutturosa"         "Rahicerus_campestris"      
##  [83] "Dorcatragus_megalotis"      "Madoqua_saltiana"          
##  [85] "Madoqua_kirkii"             "Litocranius_walleri"       
##  [87] "Eudorcas_rufifrons"         "Nanger_granti"             
##  [89] "Nanger_dama"                "Gazella_spekei"            
##  [91] "Gazella_dorcas"             "Gazella_gazella"           
##  [93] "Gazella_bennettii"          "Gazella_leptoceros"        
##  [95] "Gazella_cuvieri"            "Redunca_fulvorufula"       
##  [97] "Redunca_arundinum"          "Pelea_capreolus"           
##  [99] "Kobus_leche"                "Oreotragus_oreotragus"     
## [101] "Philantomba_maxwelli"       "Sylvicapra_grimmia"        
## [103] "Cephalophus_dorsalis"       "Cephalophus_jentinki"      
## [105] "Cephalophus_silvicultor"    "Cephalophus_spadix"        
## [107] "Cephalophus_adersi"         "Cephalophus_ogilbyi"       
## [109] "Cephalophus_callipygus"     "Cephalophus_leucogaster"   
## [111] "Cephalophus_natalensis"     "Damaliscus_pygargus"       
## [113] "Alcelaphus_buselaphus"      "Connochaetes_gnou"         
## [115] "Hippotragus_equinus"        "Hippotragus_niger"         
## [117] "Addax_nasomaculatus"        "Oryx_beisa"                
## [119] "Oryx_dammah"                "Pantholops_hodgsonii"      
## [121] "Capricornis_crispus"        "Capricornis_milneedwardsii"
## [123] "Capricornis_swinhoei"       "Naemorhedus_baileyi"       
## [125] "Naemorhedus_griseus"        "Ovis_aries"                
## [127] "Ammotragus_lervia"          "Arabitragus_jayakari"      
## [129] "Rupicapra_pyrenaica"        "Budorcas_taxicolor"        
## [131] "Pseudois_nayaur"            "Capra_sibirica"            
## [133] "Capra_nubiana"              "Capra_pyrenaica"           
## [135] "Capra_ibex"                 "Capra_caucasica"           
## [137] "Capra_falconeri"
#remove from each dataset after removing incompletes in dataset
m <- m[complete.cases(m),]
m.tree <- m[m$taxon %in% newNameList,]
df.tree <- df[rownames(df) %in% newNameList,]
testTree_crop <- drop.tip(testTree_crop,tip = testTree_crop$tip.label[!testTree_crop$tip.label %in% m.tree$taxon])

artTree <- testTree_crop
#align so that tips and m.tree are in same order
m.tree <- m.tree[artTree$tip.label,]
df.tree <- df.tree[artTree$tip.label,]
plot(artTree, cex = 0.4)

df.phylo.pca <- phyl.pca(artTree, df.tree, method="BM", mode="cov")
summary(df.phylo.pca)
## Importance of components:
##                              PC1        PC2         PC3         PC4
## Standard deviation     6.0083652 0.63086636 0.389308403 0.251244557
## Proportion of Variance 0.9816734 0.01082254 0.004121373 0.001716515
## Cumulative Proportion  0.9816734 0.99249599 0.996617358 0.998333874
##                                 PC5          PC6         PC7
## Standard deviation     0.1896223710 0.1385120336 0.078285155
## Proportion of Variance 0.0009777629 0.0005217103 0.000166653
## Cumulative Proportion  0.9993116367 0.9998333470 1.000000000
plot(df.phylo.pca)

#plot(df.phylo.pca$S[,1],df.phylo.pca$S[,2], col = habitat.colors[m.tree$habitat], 
 #    pch =16)
#text(df.phylo.pca$S[,1],df.phylo.pca$S[,2], labels=artTree$tip.label, cex = 0.5)
require(maptools)
## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
plot(df.phylo.pca$S[,1],df.phylo.pca$S[,2], col = habitat.colors[m.tree$habitat], 
     pch =16, xlab = "PC1 (98.73%)", ylab = "PC2 (0.67 %)")

#points(df.phylo.pca$L[,1:2])
#pointLabel(df.phylo.pca$L[,2:3],             # set position of labels
#           labels=rownames(df.phylo.pca$L),  # print labels
#           cex=0.75 )
plot(df.phylo.pca$S[,2],df.phylo.pca$S[,3], col = habitat.colors[m.tree$habitat], 
     pch =16, xlab = "PC2 (0.67%)", ylab = "PC3 (0.45%)")

#points(df.phylo.pca$L[,2:3])
#pointLabel(df.phylo.pca$L[,2:3],             # set position of labels
#           labels=rownames(df.phylo.pca$L),  # print labels
#           cex=0.75 )
#pointLabel(df.phylo.pca$S[,2:3],             # set position of labels
 #          labels=rownames(df.phylo.pca$S),  # print labels
 #         cex=0.75)
#loadings
#quartz(height = 6, width = 10)
par(mfrow=c(2,2))

plot(df.phylo.pca$L[,1],df.phylo.pca$L[,2], xlab="PC1", ylab="PC2", main = "Phylogenetic PCA")
pointLabel(df.phylo.pca$L[,1:2],             # set position of labels
           labels=rownames(df.phylo.pca$L),  # print labels
           cex=0.75 )
plot(df.phylo.pca$L[,2],df.phylo.pca$L[,3],xlab="PC2", ylab="PC3")
pointLabel(df.phylo.pca$L[,2:3],             # set position of labels
           labels=rownames(df.phylo.pca$L),  # print labels
           cex=0.75 )